namespace MSL {


      bool PQMRBiCGStab(const DistributedDenseMatrix< double > & A, const DistributedDenseVector< double > & b,
            DistributedDenseVector< double > & x, double & residualNorm, const double & eps) {
            
            // All (I think) variables
            DistributedDenseVector< double > p(b);
            DistributedDenseVector< double > r(p);
            DistributedDenseVector< double > rtilde(r);
            DistributedDenseVector< double > v(A.np, A.n, 1, A.r, A.c, 0.0);
            DistributedDenseVector< double > t(A.np, A.n, 1, A.r, A.c, 0.0);
            DistributedDenseVector< double > q(A.np, A.n, 1, A.r, A.c, 0.0);
            DistributedDenseVector< double > s(A.np, A.n, 1, A.r, A.c, 0.0);
            DistributedDenseVector< double > u(A.np, A.n, 1, A.r, A.c, 0.0);
            double d, rho, rhoold, tau, theta, eta;
            double alpha, omega;
            double delta, deltatilde, Phi, phi;
            double phitilde, phitildeold, gamma, pi, chi;
            double thetatilde, c, tautilde, etatilde;
            double beta, ksi, dtilde, xtilde;
            
            // Initialization
            MatVec(-1.0, A, x, 1.0, b);
            r = p;
            rtilde = r;
            MatVec(1.0, A, r, 0.0, v);
            t = v;
            d = 0.0;
            rho = dot(rtilde, r);
            tau = nrm2(r);
            theta = 0.0;
            eta = theta;
            alpha = 1.0;
            omega = 1.0;
            MatVec(1.0, A, v, 0.0, q);
            residualNorm = tau;
            
            // Iterations
            while (residualNorm >= eps) {
            
                  gaxpy(alpha, v, s);
                  gaxpy(-1.0 * alpha, q, t);
                  gaxpy(-1.0 * omega, u, t);
                  MatVec(1.0, A, t, 0.0, u);
                  
                  delta = dot(s, s);
                  deltatilde = dot(s, t);
                  Phi = dot(t, t);
                  phi = dot(rtilde, t);
                  phitildeold = phitilde;
                  phitilde = dot(rtilde, q);
                  gamma = dot(rtilde, u);
                  pi = dot(rtilde, s);
                  chi = dot(rtilde, v);
                  
                  thetatilde = std::sqrt(delta) / tau;
                  c = 1.0 / std::sqrt(1.0 + (thetatilde * thetatilde));
                  tautilde = tau * thetatilde * c;
                  etatilde = c * c * alpha;
                  omega = delta / deltatilde;
                  rhoold = rho;
                  rho = pi - (omega * phi);
                  beta = (rho * alpha) / (rhoold * omega);
                  chi = rho - (omega * gamma) + (beta * (chi - (omega * phitildeold)));
                  ksi = delta - (2.0 * omega * delta) + (omega * omega * phi);
                  alpha = rho / chi;
                  theta = std::sqrt(ksi) / tautilde;
                  c = 1.0 / std::sqrt(1.0 + (theta * theta));
                  tau = tautilde * theta * c;
                  eta = c * c * omega;
                  gaxpy((theta * theta * eta) / alpha, d, p);
                  gaxpy(etatilde, d, x);
                  r = s;
                  gaxpy(-1.0 * omega, t, r);
                  gaxpy(-1.0 * omega, v, p);
                  scale(beta, p);
                  gaxpy(1.0, r, p);
                  gaxpy(-1.0 * omega, q, v);
                  scale(-1.0 * omega, u);
                  gaxpy(beta, v, u);
                  v = t;
                  gaxpy(1.0, u, v);
                  scale(-1.0 / omega, u);
                  MatVec(1.0, A, v, 0.0, q);
                  d = s;
                  gaxpy((thetatilde * thetatilde * etatilde) / omega, dtilde, d);
                  x = xtilde;
                  gaxpy(eta, d, x);
                  r = b;
                  MatVec(-1.0, A, x, 1.0, r);
                  residualNorm = nrm2(r);
            
            }
            
            return true;
      
      }
      
      
      IBiCGStab(const DistributedDenseMatrix< double > & A, const DistributedDenseVector< double > & b,
            DistributedDenseVector< double > & x, double & residualNorm, const double & eps) {
            
            // All (I think) variables
            DistributedDenseVector< double > r(A.np, A.n, 1, A.r, A.c, 0.0);
            DistributedDenseVector< double > u(A.np, A.n, 1, A.r, A.c, 0.0);
            DistributedDenseVector< double > f(A.np, A.n, 1, A.r, A.c, 0.0);
            DistributedDenseVector< double > q(A.np, A.n, 1, A.r, A.c, 0.0);
            DistributedDenseVector< double > v(A.np, A.n, 1, A.r, A.c, 0.0);
            DistributedDenseVector< double > s(A.np, A.n, 1, A.r, A.c, 0.0);
            DistributedDenseVector< double > t(A.np, A.n, 1, A.r, A.c, 0.0);
            DistributedDenseVector< double > x(A.np, A.n, 1, A.r, A.c, 0.0);
            DistributedDenseVector< double > z(A.np, A.n, 1, A.r, A.c, 0.0);
            double sigma = 0.0;
            double pi = 0.0;
            double theta = 0.0;
            double tau = 0.0;
            double sigmaold = 0.0, sigma;
            double rhoold = 0.0, rho;
            double alpha = 0.0;
            double omega = 0.0;
            double delta = 0.0, beta = 0.0, kappa = 0.0;

            // Initialize
            r = b;
            MatVec(-1.0, A, x, 1.0, r);
            MatVec(1.0, A, r, 0.0, u);
            MatVec(1.0, *(A.getTranspose()), r, 0.0, f);
            sigma = dot(r, u);
            residualNorm = nrm2(r);
            
            // Iterations
            while (residualNorm >= eps) {
            
                  rho = phi - (omega * sigmaold) + (omega * alpha * pi);
                  delta = rho * alpha / rhoold;
                  beta = delta / omega;
                  tau = omega + (beta * tau) - (delta * pi);
                  alpha = rho / tau;
                  v = u;
                  gaxpy(beta, v, u);
                  gaxpy(-1.0 * delta, q, v);
                  MatVec(1.0, A, v, 0.0, q);
                  s = r;
                  gaxpy(-1.0 * alpha, v, s);
                  t = u;
                  gaxpy(-1.0 * alpha, q, t);
                  z = r;
                  scale(beta, z);
                  gaxpy(alpha, r, z);
                  gaxpy(-1.0 * alpha * delta, v, z);
                  phi = dot(r, s);
                  pi = dot(r, q);
                  gamma = dot(f, s);
                  eta = dot(f, t);
                  theta = dot(s, t);
                  kappa = nrm2(t);
                  omega = delta / kappa;
                  sigma = gamma - (omega * eta);
                  r = s;
                  gaxpy(-1.0 * omega, t, r);
                  gaxpy(1.0, z, x);
                  gaxpy(omega, s, x);
                  r = b;
                  MatVec(-1.0, A, x, 1.0, r);
                  residualNorm = nrm2(r);
                  MatVec(1.0, A, r, 0.0, u);
            
            }
            
            return true;
            
      }
      
      
}

